Renormalized Second-order Perturbation Theory for the Electron Correlation Energy: 

Concept, Implementation, and Benchmarks 
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We present a renormalized second-order perturbation tlieory (rPT2), based on a Kohn-Sliam (KS) 
reference state, for the electron correlation energy that includes the random-phase approximation 
(RPA), second-order screened exchange (SOSEX), and renormalized single excitations (rSE). These 
three terms all involve a summation of certain types of diagrams to infinite order, and can be viewed 
as '^renormalization" of the 2nd-order direct, exchange, and single excitation (SE) terms of Rayleigh- 
Schrodinger perturbation theory based on an KS reference. In this work we establish the concept 
of rPT2 and present the numerical details of our SOSEX and rSE implementations. A preliminary 
version of rPT2, in which the renormalized SE (rSE) contribution was treated approximately, has 
already been benchmarked for molecular atomization energies and chemical reaction barrier heights 
and shows a well balanced performance [Paier et al, New J. Phys. 14, 043002 (2012)]. In this work, 
we present a refined version of rPT2, in which we evaluate the rSE series of diagrams rigorously. 
We then extend the benchmark studies to non-covalent interactions, including the rare-gas dimers, 
and the S22 and S66 test sets. Despite some remaining shortcomings, we conclude that rPT2 gives 
an overall satisfactory performance across different chemical environments, and is a promising step 
towards a generally applicable electronic structure approach. 
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I. INTRODUCTION 

Density-functional theory (DFT)^'^ has played a sig- 
nificant role in first-principles electronic-structure calcu- 
lations in physics, chemistry, materials science, and bio- 
physics over the past decades. DFT offers an in prin- 
ciple exact formalism for computing ground-state ener- 
gies of electronic systems, but in practice the exchange- 
correlation (XC) energy functional has to be approxi- 
mated. Existing approximations to the XC functional 
can be classified into different rungs according to a 
hierarchical scheme known as "Jacob's ladder" The 
random-phase approximation (RPA))^i^ which in the con- 
text of DFT— amounts to treating the exchange en- 
ergy exactly and the correlation energy at the level of 
RPA, is on the fifth and highest rung of this ladder. 
RPA has received considerable attention (for two re- 
cent reviews, see Refs. H and |^ since its first applica- 
tion to realistic systemsJ^ This is largely due to the 
fact that RPA has shown great promise in resolving 
difficulties encountered by the local-density and gener- 
alized gradient approximations (LDA/GGAs) to DFT. 
The resolution of the "CO adsorption puzzle" jii"— the 
encouraging behavior for the "strongly correlated" /- 
electron Ce metal^i^ and the excellent performance of 
RPA (and its variants) across a wide range of systems 
including solids j^i^ii^ van der Waals (vdW) bonded 
molecules fi2r— and thermochemistry^ are just a few ex- 
amples. 

Quantitatively, however, RPA itself docs not always 
provide the desired accuracy. It was found empirically 
that the common practice of evaluating both the exact- 
exchange and the RPA correlation energy in a post- 
processing way using Kohn-Sham (KS) or generalized 
KS orbitals leads to a systematic underestimation of 
bond strengths in both molecules and solids i^°i^^'^°i^^ It- 



erating RPA to self-consistency does not alleviate this 
problem)^ Various attempts have been made in the 
past to improve the standard RPA schemej ^^i^^'^^i^^'^^" — 
with varying degrees of success. Here we will focus 
on two flavors of beyond-RPA schemes that both al- 
leviate the underbinding problem of RPA: the second- 
order screened exchange (SOSEX) ^^i^^i'^^ and the single- 
excitation (SE) correction^ SOSEX was originally for- 
mulated in the context of coupled cluster theory,— 
and accounts for the antisymmetric nature of the many- 
electron wave function. Like RPA, it can be interpreted 
as an infinite summation of a set of topologically similar 
diagrams^i^^ Adding SOSEX to RPA makes the the- 
ory one-electron "self-correlation" free. The SE correc- 
tion, on the other hand, accounts for the fact that the KS 
orbitals are not optimal for a post-processing perturba- 
tion treatment at the exact-exchange leveli^ In analogy 
to RPA and SOSEX, one can also identify a sequence 
of single excitation diagrams. Summing these to infi- 
nite order yields what we called the renormalized single- 
excitation (rSE) contribution^^ to the electron correla- 
tion energy. Combining all three contributions - RPA, 
SOSEX, and rSE - leads to the "RPA-f-SOSEX-hrSE" 
scheme, or as we shall refer to it in this work: renor- 
malized 2nd- order perturbation theory, in short rPT2 
(note that in Ref. Q we used the acronym r2PT). The 
name is inspired by second-order Rayleigh-Schrodinger 
perturbation theory (RSPT) that becomes renormalized 
through the infinite summations. This can be compared 
to the commonly used second-order M0ller-Plesset (MP2) 
method, which is a straight (bare) second-order RSPT 
based on the Hartree-Fock reference. 

A preliminary version of rPT2, in which an approx- 
imate treatment of rSE was invoked, had been bench- 
marked for atomization energies of molecules and chem- 
ical reaction barrier heights in Ref. H^- We found that 
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rPT2 gives the "most balanced" performance compared 
to other RPA-based schemes. However, this approximate 
treatment of rSE turns out to be problematic for weak 
interactions and exhibits an unphysical behavior in, e.g., 
the binding energy curve of rare gas dimers. In this work, 
we will show how a rigorous evaluation of rSE can be 
carried out. From here on, rPT2 will refer to this re- 
vised scheme and not the approximate version presented 
in Ref. H^- We will, in particular, examine the perfor- 
mance of rPT2 for weakly bonded molecules, including 
rare-gas dimers, and the widely used S22 and S66 test 
sets of Hobza and co-authorsi^"— For completeness, we 
will also revisit the benchmark sets for the G2 atomiza- 
tion energies of Curtiss et aZ.— and the chemical reaction 
barrier heights of Truhlar and co-authors22ii2 for which 
the performance of the preliminary rPT2 version was first 
tested in Ref. [si]. In addition to the concept of rPT2 
and benchmark studies, we will also present a different 
way of formulating the SOSEX term, that corresponds 
to the adiabatic connection formulation of SOSEX (AC- 
SOSEX) by Jansen, Liu, and Anygan (JLA)^ and that 
reflects our actual implementation. Our benchmark stud- 
ies show that rPT2 represents an overall improvement 
over RPA, and gives a gratifying performance across dif- 
ferent electronic and chemical environments. We also 
identify remaining shortcomings that will guide further 
developments of the theory. 

The remainder of the paper is organized as follows: In 
Sec. ini the basic theory and implementation of rPT2 is 
presented. This is followed by a systematic benchmark 
test for rPT2 for a range of systems in Sec. IIIII Conclu- 
sions are drawn in Sec. IIVI Further details of our imple- 
mentation and derivations will be given in Appendices. 



II. THEORY 

In this section the theoretical foundation of rPT2 will 
be presented. We first recapitulate the basics of the 
RPA+SOSEX method in Sec. Ill A[ and present the the- 
ory in a way that reflects its implementation in the Fritz 
Haber Institute ab initio molecular simulations (FHI- 
aims) code package . ^^i^'^ This is followed by the deriva- 
tion of an algebraic expression for the rSE term - the 
third ingredient in rPT2. A discussion of the underlying 
physics behind the rPT2 method is then presented from 
a diagrammatic point of view in Sec. Ill CI 



A. The RPA+SOSEX method 

The RPA method can be formulated in different ways 
(for a review, see Ref.iandi). In the DFT context, RPA 
can be derived from the adiabatic-connection fluctuation- 
dissipation (ACFD) thcorem^^ii whereby the RPA corre- 



lation energy is expressed as 

Y J ^^J JJ drdr'v{r,r')x 
[xT^{r\r,tuj)-xo{r',r,tu;)] . (1) 



E. 



RPA 



Xo(*w) is the KS independent-particle density-response 
function 



Xo(r,r',iw) = 



V*(r)Va(r)V^(r')V'a(r') 



-I- c.c. 



■ ea-iui 



(2) 



where '!/'i.a(r) and e^^a are the KS single-particle orbitals 
and orbital energies, and c.c. the "complex conjugate". 
Here and in the following we adopt the following conven- 
tion: I, j correspond to occupied and a, & to unoccupied 
(or virtual) spin orbitals, whereas p, q apply to general 
cases. Xa^^I**^) ^1- © RPA response function 

of a fictitious system with a scaled Coulomb interaction 
l^^pj (with < A < 1), and satisfies the Dyson equation 

XA = Xo + XoAwxA ■ (3) 

Representing xo and v in the "particle-hole basis" 
{i/'*(r)i/)a(r), V'a(r)V'i(i')}j one can obtain the RPA cor- 
relation energy by solving the following eigenvalue 
probleroiS 



A 

-B 



B* 

-A* 



(4) 



where Aiajb = {e-a - ^i)6ij6ab + {ib\aj), and Bia.jb 
{ij\ab). The two-electron Coulomb integrals are 



{pq\rs) 



%{xi)ll)r{xi)i)*{x2)'lps{x2) 

dxidx2 j ] , (5) 



ri - r2 



where x ~ (r, a) is a combined space-spin variable. As 
demonstrated by Furche^ after solving Eq. ([4]), the RPA 
correlation energy can be written as 



E. 



RPA 



1 



Tr(a; - A) 



^ w„ - ^ Ai, 



, (6) 



where implies that the summation over n is restricted 
to positive eigenvalues a;„. 

Scuseria et al. demonstrated that an equivalent for- 
mulation of the RPA correlation energy of Eq. ^ can 
be obtained from an approximate coupled-cluster dou- 
bles (CCD) theory^ in which only the "ring diagrams" 
are kept (see the first row of Fig.[T|). In the CCD theory, 
only double excitation contributions are included in the 
"cluster operator" which generates the interacting many- 
body ground-state wavefunction through the exponential 
ansatz. By contrast, in the more often used CCSD ap- 
proach, both single and double excitations are included. 
Within the CCD formulation of RPA, the key quantities 
are the {direct) ring-CCD amplitudes Tiajb, which (in 





FIG. 1: Goldstone diagrams for RPA (first row) and SOSEX (second row) contributions. Dashed lines represent bare Coulomb 
interactions, and full lines correspond to KS electrons (arrow up) and holes (arrow down). Third row: RPA+SOSEX energy 
in the coupled-cluster context. The wiggly together with the arrowed, solid lines represent the direct ring-CCD amplitudes 
Tiajb (see Eq. ((9)1). The contraction between the direct ring-CCD amplitudes and the bare Coulomb interactions (dashed lines) 
yields the RPA+SOSEX correlation energy. 



the case of real canonical spin orbitals) are determined 
by the following Riecati equation, 



B + AT + TA + TBT = . 



(7) 



Due to the quadratic nature of this equation, one should 
take care to ensure that the physical solution is taken4^ 
The RPA correlation energy in this ring-CCD formula- 
tion is then given by 



EY''^\Tr{BT) = \Y.^^3\ab)T, 



jb^ia ■ 



(8) 



We note that this is often called direct RPA in the 
quantum chemistry literature to emphasize the fact that 
higher-order exchange-type contributions arc not in- 
cluded. 

Now the RPA+SOSEX correlation energy can be 
conveniently introduced^^'^'^ by antisymmetrizing the 
Coulomb integral in Eq. 



^RPA+SOSEX ^lTr(Br) = i^(^,||a6)T, 



jb,ia ■ 



(9) 



where B.,a,jb {ij\\ab) ^ {ij\ab) - {ij\ba). The SOSEX 
correction term itself is 



E^^osE^^iBT) = -lY.(^j\ba)T,,,^. 



(10) 



Physically, the SOSEX correction introduces higher- 
order exchange processes that can also be represented by 
an infinite summation of Goldstone diagrams (see the sec- 
ond row of Fig. H]). This infinite summation is condensed 
into the ring-CCD amplitudes whose contraction with 
the bare Coulomb interaction (after antisymmetrization) 
yields the RPA+SOSEX correlation energy as illustrated 
by the third-row diagrams in Fig. [T] 

In a coupled cluster code the SOSEX energy can be 
readily computed once the direct ring-CCD amplitudes 
Tia.jb are available. A slightly different variant of SO- 
SEX can be obtained in the ACFD framework, as shown 
by JLAii. We will show later, that although not identi- 
cal, these two SOSEX formulations produce very similar 
results. Our implementation in the FHI-aims code^i^ 
follows the ACFD route. To illustrate our approach let 
us first present an alternative way to Eq. ^ of express- 
ing the RPA correlation energy within ACFD before we 
introduce the corresponding SOSEX extension. Eq. ^ 
yields 



(iw) =Xo(*w) + Xo(«'^)AwXa^^(«'^) 
=Xo(«w) + Axo(«^^)i'Xo(«w)+ 



(11) 



The RPA correlation energy in Eq. ([T]) can then be rewrit- 
ten as 



ij,ab 



1 



E^^^ = -— J d.^ I '^'^Tr [xo(«w)i;xo(«w) • Aw + xo(iw)uxo('iw) ■ A^wXo(»t^)f H J (12) 







= y" dA y" du}Tr[xoiiu})vxoiii^)Wx{iu})] 

= -2^y duTr[xo{iuj)vxo{t^)W{iuj)] , (13) 
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where 



Wxiiuj) = Aw/(1 - Xxo{iuj)v) 



(14) 



is the coupUng-constant-dependent screened Coulomb in- 
teraction and 



W{iu}) = / dXWx{iu}) 
Jo 



(15) 



the couphng-constant-averaged screened Coulomb inter- 



action. In this context we would like to point out that 
the first diagram in the third row of Fig. [T] can alter- 
natively be interpreted as the pictorial representation of 
equation Now the bubbles correspond to xo, dashed 
lines to the bare Coulomb interaction, and wiggly lines 
to the corresponding screened interaction W{iuj). 

Expressing xo again in terms of the "particle-hole ba- 
sis" (defined below Eq. ©) and using Eq Eq. ^ 
can be recast into 



J 



ia.jb 



{aj\ib){ib\W{iuj)\aj) 



{ab\ij){ij\W{ioj)\ab) 



(ej - Ca - iuj)iej - £& - »w) (cj - ea - iw)(ej - et + i^^) 

{ij\ab){ab\Wiiuj)\ij) ^ {ib\aj){aj\W{iLu)\ib) 
(e,; - Eq + iuj){ej - Efc - iuj) (e., - Ca + iw)(ej - £6 + iw) 



(16) 



where {pq\W{iLL:)\rs) is defined in analogy to {pq\\7\s) in 
Eq. ([5]), by replacing the bare Coulomb interaction v by 
the screened (and frequency-dependent) one, W{iuj). 

For real canonical spin orbitals we find {aj\ib) ~ 
{ib\aj) = {ab\ij) = {ij\ab). The same relations hold 
for the screened Coulomb repulsion integrals. The above 
equation then simplifies to 

1 l"^ 

Fia{iuj)F-jb{ii^) (17) 

with the factors 

F^a{ioj) = 2(e, - ea)/[(e, - Eaf + uj^] . (18) 

Now, in analogy to the [direct) ring-CCD formulation 
of SOSEX in Eq. (fTO|l. one can obtain a corresponding 
SOSEX term (the so-called "AC-SOSEX") from Eq. (fTT)) . 
by exchanging the "a, 6" indices in {ij\ba) (with an addi- 
tional minus sign), 

^AC-SOSEX ^ _ / dujy2{ij\ba){ij\W{iLo)\ab)x 

Fia{ii^)T.jb{iuj) . (19) 

Then, using the resolution-of-identity technique )"^^ ,47- _49 
Eq. p9p can be implemented with relative ease. The 
implementation details of Eq. p9|) in FHI-aims are pre- 
sented in Appendix El 



To make closer contact with the expression given in 
Ref. Hll, we note that Eq. (|17p can be further rewritten: 

i^c^"-™ = -^E(*^>«)^-.'" (20) 

where 

1 /■°° 

Pia,ib = - duj{ij\W{iuj)\ab)Tia{iuj)J^jb{tuj) (21) 

Jo 

is the coupling-strength averaged (two-particle) density 
matrix. 

As shown by JLAM^, Eq. ([20l) is usually not identical 
to the original ring-CCD based SOSEX in Eq. ([TOl) (ex- 
cept for one- and two-electron cases). However, the dif- 
ference between them is very small (relative difference in 
RPA-f SOSEX correlation energy less that 0.15%), as first 
noted in Ref. [s^ and also confirmed here. In table H] we 
present the RPA and SOSEX correlation energies {Ef-^^ 
and E^'^^^^), as weh as the RPA and RPA-I-SOSEX at- 
omization energies for five molecules. The vanishingly 
small differences in the RPA energies are due to the dif- 
ferent implementations in FHI-aims and the development 
version of the GAUSSIAN^ code (e.g., FHI-aims em- 
ploys the RI approximation and treats the Gaussian or- 
bitals numerically). The difference in the SOSEX and 
AC-SOSEX correlation energies reflects the intrinsic dif- 
ferences of the two SOSEX formulations. Nevertheless, 
the differences are very small and have little practical 
importance, in particular for atomization energies. 



Our benchmark results presented in section IIIII are 



based on the AC-SOSEX scheme. However, since the 
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TABLE I: RPA and SOSEX (total) correlation energies (in Hartree), as well as RPA and RPA+SOSEX atomization energies 
(in kcal/mol) for five molecules. The "AC-SOSEX" numbers are computed using FHI-aims based on Eq. (|19p . whereas the 
original ring-CCD based SOSEX numbers are computed using a development version of the GAUSSIAN^° suite of programs. 
All calculations were done with Gaussian cc-pVQZ basis set and frozen-core (Is) approximation. The reference orbitals are 
obtained using the CCA functional constructed by Perdew, Burke, and Ernzerhof (PBE).— Note that in the upper part of the 
table only the RPA or (AC-)SOSEX correlation contribution is included, whereas in the lower part the numbers are obtained 
from the total energy (including also the Hartree-Fock part) differences. 









Correlation enerj 


Ij (Hartree) 










RPA 






AC-SOSEX/SOSEX 






FHI-aims 


GAUSSIAN 


difference 


FHI-aims 


GAUSSIAN 


difference 










(AC-SOSEX) 


(SOSEX) 




CO 


-0.593778 


-0.593786 


0.000008 


0.218954 


0.217977 


0.000977 


N2 


-0.606368 


-0.606391 


0.000023 


0.224069 


0.222955 


0.001114 


02 


-0.730348 


-0.730364 


0.000016 


0.283384 


0.281073 


0.002311 


CH4 


-0.381735 


-0.381730 


-0.000005 


0.155242 


0.154933 


0.000309 


C2H2 


-0.539435 


-0.539439 


0.000006 


0.207348 


0.206514 


0.000834 








Atomization enerf 


;y (kcal/mol) 










RPA 




RPA+ 


AC-SOSEX/RPA-KSOSEX 




FHI-aims 


GAUSSIAN 


difference 


FHI-aims 


GAUSSIAN 


difference 










(AC-SOSEX) 


(SOSEX) 




CO 


239.16 


239.18 


-0.02 


246.88 


246.86 


0.02 


N2 


217.58 


217.59 


-0.01 


209.24 


209.10 


0.14 


02 


108.02 


108.03 


-0.01 


98.11 


98.71 


-0.60 


CH4 


400.15 


400.13 


0.02 


415.33 


415.29 


0.04 


C2H2 


373.43 


373.45 


-0.02 


391.63 


391.71 


-0.08 



numerical difference between the two SOSEX flavors are 
very small, our conclusion should also apply to the orig- 
inal ring-CCD based SOSEX. 



B. The rSE correction and the 
semi-canonicalization method 

In Rcf. [23, we showed that it is advantageous to 
complement the RPA correlation energy with a cor- 
rection term arising from single excitations. The sin- 
gle excitation correction derives directly from Rayleigh- 
Schrodinger perturbation theory (RSPT) and adopts a 
simple form in terms of the single-particle orbitals 



E: 



SE 



E 



E 



\fia I 



(22) 



Here |'0i(a)) ^-nd ^i(a) refer to occupied (unoccupied) 
Kohn-Sham (KS) orbitals and the corresponding orbital 
energies. / is the single-particle Hartree-Fock (HF) 
Hamiltonian, or the so-called Fock operator. We have 
presented the derivation of Eq. (|22|) already in Ref. [20l . 
but include it here for completeness in Appendix [B) De- 
noting the single-particle KS Hamiltonian , we ob- 
tain {'^p^\f\i)a) = {'^p^\h^ + M\l}:a) = {1p^\l^v\^Pa) whcu 

ipi^ ipa are eigcnfunctions of hP . Au is the difference 



between the HF exact-exchange potential and the KS 
exchange-correlation potential. A similar SE contribu- 
tion is encountered in the context of KS density func- 
tional perturbation theoryi^"— However, we emphasize 
that here we followed the procedure of RSPT to derive 
Eq. p2p . instead of the ACFD formalism, which requires 
the electron-density to be fixed along the adiabatic- 
connection path. Whether the two procedures will yield 
significantly different results is a subject of further stud- 



From the viewpoint of RSPT, Eq. ([22]) represents a 
second-order correlation energy. As such it suffers from 
the same divergence problem as 2nd-order M0ller-Plesset 
perturbation theory for metallic systems when the single- 
particle KS gap closes. A remedy suggested in Ref. [20 
was to follow the RPA spirit and to sum a sequence of 
higher-order SE terms to infinite order. Such higher- 
order SE terms can also be represented in terms of Gold- 
stone diagrams, as illustrated in Fig. [2l We refer to this 
infinite summation of SE terms as renormalized single 
excitations (rSE) as alluded to in the introduction. 

The influence of the rSE correction was first examined 
in Ref. Is^ . albeit in an approximate way. There a so- 
called "diagonal" approximation to rSE (denoted here as 
"rSE-diag") was used, in which only terms with "i — j — 
k = • • ■ " and "a = 6 = c = • • • " were included. The 
remaining "off-diagonal" terms were omitted. A similar 
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4th-order 



FIG. 2: Goldstone diagrams for a sequence of correlation-energy terms arising from single excitations. Summing these up to 
infinite order yields the renormalized single excitation (rSE) contribution. Here Avpq — {tpp\f — and note Avia = fia- 



approximation has been used in summing up the Epstein- 
Nesbet ladder-type diagrams in Ref. [SJ. In this way, 
the sequence of diagrams falls into a geometrical scries. 
Summing them up yields the following simple expression 



one has 



E: 



rSE-diag 



fia I " 



Avii - AVa 



(23) 



where Avpq = (tpplAvlipq). The additional term Avu — 
Avaa that appears in the denominator is negative defi- 
nite and removes the divergence problem even for van- 
ishing KS gaps. The addition of rSE-diag to RPA and 
RPA+SOSEX has been bcnchmarked for atomization en- 
ergies and reaction barriers in Ref. H^. We found that 
the renormalization (i.e., going from SE to rSE-diag) has 
a tendency to slightly reduce atomization energies, but 
the overall effect is not significant. For chemical reaction 
barrier heights, on the other hand, the renormalization 
is crucial for the transition states, that typically have a 
rather small energy gap. 

The diagonal approximation in Eq. (|23|) is not invariant 
under unitary transformations in the space of occupied 
and unoccupied orbitals. More importantly, however, 
it can lead to an unphysical behavior in the potential- 
energy surface of weakly interacting systems, as will be 
shown in Sec. IIII A II Recently we discovered that it is 
straightforward to include the "off-diagonal" elements as 
well, and to treat the rSE term rigorously. In Appendix [Cl 
we illustrate in detail how the infinite summation of the 
diagrams depicted in Fig. [2] can be carried out. Here we 
only present the key steps that lead to the final expres- 
sion, and that arc needed in practical calculations. 

First, the occupied and unoccupied blocks of the Fock 
matrix (evaluated with KS orbitals) need to be con- 
structed 

fij = {ipi\f\i>j) = ^tSij + Avij 

fab = (V'al/IV'fc) <^aSab + AVab ■ 

The second step is to diagonalize the fij and the fab block 
separately. Denoting the eigenvector matrices as O and 



J2f^kOkj^O.,ej 

k 

faMcb = Uabh , 



(24) 



where ej and eb are the eigenvalues of the occupied and 
unoccupied blocks of the Fock matrix, respectively. This 
procedure is known as semi-canonicalization in quantum 
chemistry (see e.g. Ref. Issl ). The final rSE expression, 
equivalent to the infinite-order diagrammatic summation, 
is given by 



E. 



rSE 



\fl 



(25) 



where fia correspond to the "transformed" off-diagonal 
block of the Fock matrix 



ijl^*abfjb ■ 

jb 



(26) 



This is a surprisingly simple result: the final rSE expres- 
sion is formally identical to the 2nd-order SE one; only 
that the meaning of the energy eigenvalues and the tran- 
sition amplitudes has to be modified. The equivalence 
of Eq. ((25|) to the algebraical expression from a direct 
evaluation of the diagrams in Fig. [2] is demonstrated in 
appendix [Cl 



C. The concept of rPT2 viewed from its 
diagrammatic representation 

Initially the RPA+SOSEX and RPA-|-(r)SE schemes 
were developed separately22i2£ in an effort to improve the 
accuracy of the RPA method. In Ref. [13 it was found that 
adding both terms to RPA leads to even better accuracy 
in general, and that the combined RPA-|-SOSEX-|-rSE 
(=rPT2) scheme represents the most balanced approach 
for describing both atomization energies and reaction 
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2nd-order 3rd-order 



FIG. 3: (Color online) rPT2 represented in terms of Goldstone 
diagrams. The three rows of (infinitely summed) diagrams 
represent the three components of rPT2: RPA, SOSEX, and 
rSE. The first column shows the (only) three terms in nor- 
mal (bare) 2nd-order Rayleigh-Schrodinger perturbation the- 
ory based on a KS reference. 



barrier heights. To elucidate the nature of rPT2, the 
Goldstone diagrams for the three ingredients of this the- 
ory are shown together in Fig. [31 All three pieces are 
characterized by an infinite summation of diagrams with 
the same topological structure. The leading terms in the 
three series are the second-order direct (Coulomb), the 
second-order exchange, and the SE term, respectively. In 
other words, these leading terms are exactly the (only) 
three terms that one would encounter in second-order 
Rayleigh-Schrodinger perturbation theory, based on an 
(approximate) KS reference Hamiltonian. Only the SE 
term would vanish if the perturbation series were to be 
build on the HF reference. In essence, the theory is ex- 
act at second order, and for higher-order contributions we 
follow the strategy of "selective summation to infinite or- 
der" , following the spirit of the RPA. This "infinite-order 
summation" effectively renormalizes the three terms of 
the (bare) second-order perturbation theory (PT2), rep- 
resented by the the blue diagrams in Fig. [S] We expect 
the renormalized method, i.e. rPT2, to be more gen- 
erally applicable than the bare PT2, which, e.g., suffers 
from notorious divergence problems for systems with zero 
direct gapi^iS 

As a perturbation theory, rPT2 will necessarily depend 
on the reference Hamiltonian, or equivalently a set of in- 
put single-particle orbitals. In practice, rPT2 works best 
when based on Kohn-Sham Hamiltonians, that yield a 
smaller gap than generalized KS or HF ones. This is 
directly related to the fact that the underbinding er- 
ror of RPA will be even more pronounced for HF or 
generalized KS reference Hamiltonians, as evidenced by 
the significant RPA@HF error for the G2 atomization 
energies;^ and the severely underestimated RPAQHF 
(40%) Cg coefficients^ (here and in the following, we use 



"mei/iorf@re/erence" to denote which method is based on 
which reference state) . For a variety of KS Hamiltonians 
(i.e. with local, multiplicative potentials), RPA results 
were found to be insensitive to the actual choice of the 
reference Hamiltoniauj ^^i^^ In this work, we will therefore 
choose the most popular non-empirical GGA functional 
PBE as the reference; also to be consistent with our pre- 
vious workj^i2£i2i The insensitivity of RPA to reference 
KS Hamiltonians carries over to rPT2. 



III. RESULTS 

In this section we will benchmark the performance 
of rPT2 for weak interaction energies (rare-gas dimers, 
S22 and S66 test sets by Hobza and coworkers^^*^) , at- 
omization energies (from the G2-I test set by Gurtiss 
et aZi^i^), and chemical reaction barrier heights (38 
hydrogen-transfer and 38 non-hydrogen-transfer chemi- 
cal reactions of Truhlar and coworkers^SiiS) . AH cal- 
culations were performed with the local-orbital based 
all-electron FHI-aims code42ii^ As mentioned in sec- 
tion III Al the SOSEX term in this work corresponds to 
"AC-SOSEX" based on Eq. For brevity we will 

simply refer to it as SOSEX in the following. For the 
frequency integration in our RPA and SOSEX calcula- 
tions, we use a modified Gauss-Legendre grid^ with 40 
points. For the A integration in Eq. psp . we use a nor- 
mal Gauss-Legendre grid with 5 points. These settings 
guarantee sufficient accuracy for the benchmark studies 
presented in this work. The basis sets employed in the 
calculations will be specified later when discussing the 
results. Convergence tests are shown in Appendix [Pl 

A. Weak interactions 

One prominent feature of RPA-based approaches is 
that the ubiquitous vdW interactions can be captured in 
a seamless manner i^i*^ The long-range behavior of the 
RPA interaction energy between two closed-shell molec- 
ular systems decays as Ce / where the Cg value is dic- 
tated by the RPA polarizability of the monomeri^iS 
Many-body terms that go beyond the pair-wise summa- 
tion are also automatically contained in this approachi^ 

Benchmarking the performance of RPA and related 
methods for vdW bonded systems has been a very active 
enterprise] ^^i^'''^^^^^'^^!^^" — It has been demonstrated 
that the standard RPA approach exhibits a systematic 
underbinding behavior for molecules, in particular vdW 
bonded onesj^ We have previously shown that SE-type 
corrections ameliorate this problem;^ but the influence 
of the SOSEX correction has not been systematically 
benchmarked for vdW systems yet, with the exception of 
He2 and Ne2.— It is therefore interesting and timely to 
examine how rPT2, that combines both types of correc- 
tions, performs for noncovalent interactions. Some rPT2 
results for Ar2 and S22 have featured in our recent review 



on RPAi^ Here we extend the benchmark study to other 
rare-gas dimers and also the larger S66 test set. 



1. Rare-gas dimers 

First, we demonstrate the pathological behavior of 
rSE-diag for weak interactions, highlighting the impor- 
tance of including the "off-diagonal" terms in the rSE 
summation to make the theory invariant with respect to 
orbital rotations. In Fig. S] the binding energy of Ar2 is 
plotted for PBE, RPA, and RPA plus different versions of 
single excitation corrections (RPA+SE, RPA+rSE-diag, 
RPA-t-rSE). While PBE, RPA, and RPA+SE all show 
their characteristic behaviors, the behavior of RPA+rSE- 
diag is weird. The binding energy curve develops un- 
physical undulations away from equilibrium. Moreover, 
the asymptotic limit does not follow the correct 
behavior, and the curve even reaches above the energy 
zero at large bonding distances (see the inset of Fig. 
Naturally, this problem also carries over to rPT2-diag 
(not shown). It is reassuring, however, to observe that 
this pathological behavior disappears in the upgraded 
RPA+rSE scheme, which yields a binding energy curve 
in close agreement with the Tang-Toennies reference 
curve^^ obtained from a simple analytical model with ex- 
perimental equilibrium bond distance and binding energy 
as input parameters. This model can accurately repro- 
duce empirical data^ and agrees excellently with high- 
level quantum-chemical, e.g., CCSD(T) calculationS] ^°i^^ 
Coming back to the rSE discussion, the pathological be- 
havior is thus caused by the diagonal approximation, and 
not inherent to the rSE scheme itself. In the remainder of 
our discussion on weakly interacting systems, we there- 
fore only present results for the upgraded RPA+rSE and 
rPT2 schemes. 

The full set of binding energy curves for He2, Ne2, 
Ar2, and Krs obtained with PBE, MP2, RPA, rPT2, 
as well as the "intermediate" schemes RPA+rSE and 
RPA+SOSEX are then shown in Fig. El PBE does not 
contain long-range dispersion interactions by construc- 
tion, and therefore decays too fast at large separations. 
Around the equilibrium region, PBE vastly overbinds 
He2 and Ne2, and underbinds Ar2 and Kr2. MP2 shows 
the opposite trend, although it performs better at a quan- 
titative level. RPA systematically underbinds all dimers. 
This underbinding is most significant for IIe2 and Ne2. 
Adding the rSE correction leads to a substantial improve- 
ment for all dimers. With the largest available Dun- 
ning Gaussian basis sets2^ (aug-cc-pV6Z for He, Ne, Ar 
and aug-cc-pV5Z for Kr), RPA+rSE shows nearly perfect 
agreement with the reference curve for He2 , overshoots a 
little bit for Ne2, and slightly underbinds Ar2 and Kr2. 
The SOSEX correction, on the other hand, has very little 
effect on the binding energies of these purely dispersion- 
bonded systems. As a result, rPT2 lies almost on top of 
RPA+rSE. The overall accuracy of RPA+rSE and rPT2 
for rare-gas dimers is very satisfactory, in particular since 




Bond length (A) 

FIG. 4: (Color online) Binding energy curves for Ar2 com- 
puted with PBE and RPA-based approaches (standard RPA, 
RPA+SE, RPA+rSE, and RPA+rSE-diag based on PBE), in 
comparison with the Tang-Toennies reference curve. The re- 
sults are obtained using the Gaussian "aug-cc- pV6Z" ?? basis 
set. The basis set superposition error (BSSE) is corrected 
here and in all following calculations using the counterpoise 
correction scheme.— 

no adjustable parameters are used in these schemes. 



2. S22 and S66 test sets 

A widely used benchmark set for weak interactions are 
the S22 molecular complexes designed by Jurecka et al.,— 
for which accurate reference interaction energies obtained 
using the the CCSD(T) method are availablci^ This 
molecular test set includes the most common types of 
non-covalent interactions: hydrogen bonding, dispersion- 
dominated bonding, and those of mixed character. The 
performance of RPA and some of the RPA-related meth- 
ods have been benchmarked for this test seti ^^i^^'^'^i'^^ 
Similar to correlated quantum chemical methods, the 
quality of basis sets for RPA calculations is a significant 
issuejiSiiSii^ Using our numerical atomic orbital (NAO) 
tier 4 basis plus additional diffuse Gaussian functions 
from the aug-cc-pV5Z set (denoted as "iier 4 + a5Z- 
d"— ; see also Appendix |D|, we obtained a mean ab- 
solute error (MAE) of 0.90 kcal/mol in RPAOPBE for 
S22, fairly close to the 0.79 kcal/mol reported by Es- 
huis and Furcho^^ using Dunning's Gaussian basis sets 
extrapolated to the complete basis set (CBS) limit. In 
Appendix [D] the convergence behavior of these two types 
of basis sets is shown for the methane dimer. In this 
work we will continue to use the "iier 4 + a5Z-d" ba- 
sis set, bearing in mind that the absolute numbers could 
carry an uncertainty of 0.1 kcal/mol (4 meV), which will 
however not affect our discussion here. 

In Fig. [6] the relative errors from RPA+rSE, 
RPA+SOSEX, and rPT2 are presented for each indi- 
vidual molecule of the S22 set. Results from RPA and 
RPA+SE, as well as from PBE and MP2 are also in- 
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Bond length (A) 

FIG. 5: (Color online) Binding energy curves for rare-gas dimers computed with RPA-based approaches, in comparison with 
PBE, MP2, and the Tang-Toennies reference curves. He2, Ne2, and Ar2 results are obtained using the aug-cc-pV6Z basis set, 
and Kr2 using the aug-cc-pV5Z basis set. All RPA-type calculations are based on the PBE reference. 



eluded for comparison. PBE and MP2 are both per- 
forming well for hydrogen-bonded moleeules where the 
electrostatic interactions dominate, but PBE underbinds 
the dispersion-dominated and those of mixed-character 
significantly, while the opposite is true for MP2. RPA- 
based methods are performing much better than PBE 
and MP2 for these two types of interactions. RPA-|-rSE 
falls between RPA and RPA-I-SE, although it lies closer 
to RPA-I-SE. For hydrogen-bonded molecules, RPA+rSE 
improves over RPA-I-SE, with the latter overbinding 
these molecules noticeably. Moreover, it is interesting 
to note that RPA+SOSEX improves over RPA appre- 
ciably for hydrogen- and mixed-bonding, but much less 
so for dispersion-bonded molecules. This is consistent 
with its performance for rare-gas dimers. Now, com- 
bining rSE and SOSEX, rPT2 performs equally well or 
better for dispersion-dominated and mixed-bonding, but 
overshoots significantly for hydrogen-bonding. So far this 
is the only case we have found, for which combining rSE 
and SOSEX worsens the description. Finally we note 
that for TT-stacked systems like the benzene dimer (# 
11) RPA gives a substantial error, but neither rSE nor 
SOSEX significantly improves upon RPA. This warrants 
further attention in future studies. 

Recently the S22 test set has been extended to an even 
larger, more comprehensive and balanced test set called 
S66i^ This overcomes several shortcomings of S22, e.g. 
the strong bias towards nucleic-acid-like structures. We 
also performed benchmark calculations with RPA, rPT2, 
and related computational schemes for this test set, and 
the results are presented in Fig. [T] The overall perfor- 
mance for S66 is very similar to that observed for S22. 
In brief, RPA-f-rSE performs better (or slightly better) 
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FIG. 6: (Color online) The percentage errors for the S22 test 
set for RPA-derived computational schemes (based on PBE 
reference orbitals), in comparison to PBE and MP2. The 
CCSD(T)/CBS results of Takatani et ali^ are used as refer- 
ence. Lines are guides to the eye. 



than RPA-(-SE, which itself is a significant improvement 
over the standard RPA method. Adding SOSEX, the 
resultant rPT2 approach performs even (slightly) bet- 
ter than RPA-|-rSE for dispersion and mixed interac- 
tions. However, this is not the case for hydrogen bonds, 
where rPT2 clearly overshoots and the strength of hy- 
drogen bonds becomes overestimated. Overall, for weak 
interactions RPA-f-rSE outperforms other computational 
schemes benchmarked here, and yields a MAE of 10.1 
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FIG. 7: (Color online) MAEs (in both meV and kcal/mol) 
for the S66 test set given by RPA, rPT2 and related schemes 
(based on PBE reference orbitals), in addition to PBE and 
MP2. The CCSD(T) results of Rezac et alM at the CBS 
limit are used here as reference. 



FIG. 8: The MAEs (in both meV and kcal/mol) of the G2-I 
atomization energiea^S obtained with PBE, MP2, RPA, rPT2, 
and related methods. The Gaussian cc-pV6Z basis se1j^ was 
used in all calculations. Reference data are from Ref. \TA 



meV (or 0.23 kcal/mol). 



B. G2 atomization energies 

The atomization energy of molecules is a key quan- 
tity in thermochemistry. RPA has been tested for this 
quantity in early works ^^2*2^ where a pronounced under- 
binding behavior was observed. In a recent work, Paier 
et al^ reported a detailed study of the atomization ener- 
gies of the G2-I seli^ using RPA and its variants, includ- 
ing the rPT2-diag scheme as discussed before. To test 
the influence of the off-diagonal elements of rSE in the 
rPT2 scheme, we present in Fig. [8] the MAEs for RPA, 
rPT2-diag, rPT2, and related methods. Some of these 
results were already included in our recent review paper 
on RPAi^ In brief, the MAE for RPA is significantly re- 
duced when adding the (r)SE or SOSEX corrections. In 
this case, RPA-t-rSE yields a slightly larger MAE than 
RPA-I-SE. Combining the rSE and SOSEX corrections, 
rPT2 reduces the MAE further by a factor of two. In 
contrast to the nonbonded interactions discussed in the 
previous section, the difference between rPT2 and rPT2- 
diag is small (0.18 kcal/mol or 8 meV difference in MAE), 
validating our previous conclusions regarding the atom- 
ization energies in Ref. [U that were based on the rPT2- 
diag scheme. 

In this context we would like to warn that, despite 
the success of RPA-I-SOSEX and rPT2 for describing 
the atomization energies on average, adding SOSEX to 
RPA makes things worse (more underbinding) for cer- 
tain molecules (in particular O2 and N2), and this prob- 
lem also carries over to rPT2. A detailed investigation of 
this issue is beyond the scope of this paper, and will be 
carried out in future work. 



C. Barrier heights 

To complete our discussion, we address chemical re- 
action barrier heights. For this purpose we chose 
the HTBH38 and NHTBH38 test set of Truhlar and 
coworkers RPA-based methods were benchmarked 
in previous studies^'^'^, and we here revisit this set with 
the upgraded version of rPT2. The MAEs for our differ- 
ent schemes arc shown in Fig. [S] Standard RPA performs 
remarkably well for reaction barrier heights compared to 
all alternatives. This has been rationalized by Henderson 
and Scuseria4^ to be due to the inherent self-correlation 
error in RPA that mimics "static correlation" (i.e. the 
(near) degeneracy of two (or more) determinants) , lead- 
ing to an excellent description of the transition states due 
to partial error cancellation. Unfortunately, any attempt 
to correct RPA deteriorates its performance in this case. 
In particular, the RPA+SE method provides a bad de- 
scription of the transition states, resulting in errors that 
are even larger than in PBE. The RPA-t-SE error reduces 
when the SE term is renormalized in RPA+rSE. The er- 
rors in RPA-l-rSE and RPA-I-SOSEX tend to cancel each 
other, and by combining the two schemes, rPT2 gives a 
much more satisfactory description of the barrier heights. 
Similar to the G2-I test set, the difference between rPT2 
and rPT2-diag is small (0.33 kcal/mol for HTBH38 and 
0.25 kcal/mol for NHTBH38 in MAE) compared to the 
variation among other schemes. 

IV. CONCLUSIONS 

In summary, the rPT2 method comprises an infinite 
summation of three distinct series of diagrams: RPA, 
SOSEX, and rSE. As is obvious from its diagrammatic 
representation, rPT2 can be viewed as a renormalization 
of bare second-order perturbation theory - the latter be- 
ing the leading term of rPT2. In this work we derived an 
alternative way to express the SOSEX correlation energy, 
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FIG. 9: (Color online) The MAEs (in both meV and 
kcal/mol) of the HTBH38 and NHTBH38 test sets for barrier 
heights, obtained with PBE, MP2, RPA, rPT2, and related 
methods (based on PBE). Reference data are from Ref. |39||40| . 
Gaussian cc-pV6Z basis sets were used in all calculations. 



discussed in detail how to sum up the "ofF-diagonal" el- 
ements in rSE, which were neglected in previous works, 
and presented the concept of rPT2 from a diagrammatic 
point of view. We benchmarked the performance of rPT2 
and related approaches (RPA+rSE, RPA+SOSEX), fo- 
cusing on weakly interacting molecules. We found that 
rPT2 works well for dispersion and mixed-type interac- 
tions, but for hydrogen bonds it over-corrects the under- 
binding behavior of RPA. We also examined the influence 
of the previously neglected "off-diagonal" elements in the 
rSE correction and found that, for weak interactions, it is 
crucial to include them, whereas for atomization energy 
and reaction barrier heights, the off-diagonal elements 
only have a minor effect. We also found that the SO- 
SEX correction improves the description of electrostatic 
interactions substantially, but has very little effect on 
dispersion interactions. rSE, on the other hand, leads to 
a better description of both electrostatic and dispersion 
interactions. 

Overall rPT2 provides a conceptually appealing, and 
diagrammatically systematic way for going beyond RPA. 
Although it does not always deliver the best accuracy 
in every single case compared to other RPA-based ap- 
proaches, it provides the most "balanced" description 
across various different electronic and chemical environ- 
ments. We thus consider the rPT2 scheme as a natural 
step for extending and improving the RPA method. The 
successes and shortcomings of rPT2 documented in this 
work provide a useful basis for developing more accu- 
rate, robust, and generally applicable electronic structure 
methods in the coming years. 
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Appendix A: Implementation of AC-SOSEX in 
FHI-aims 

The RPA implementation in FHI-aims has been de- 
scribed in detail in Ref. Here we will give a brief 
account of the SOSEX implementation in our code. The 
energy expression that we would like to evaluate is 

£,AC-sosEX = / du^Y, {^j\ba){ij\W{iu:)\ab)x 

^ la.jb 

Tia{iuj)Tjb{iuj) (Al) 

where {ij\ba) are the two-electron Coulomb integrals de- 
fined in Eq. (O, and {ij\W{ioj)\ah) are the correspond- 
ing (coupling-constant-averaged) screened Coulomb inte- 
grals. The frequency-dependent factor J-ia{iuj) is defined 
in Eq. (fT5)) . 

In analogy to the RPA case, the basic technique to 
evaluate the two-electron integrals in our code is the 
resolution-of-identity. We chose the Coulomb metric, de- 
noted "RI-V" in the following. Here we would like to 
emphasize that "RI-V" is a highly accurate method, and 
the error incurred thereby is vanishingly small for prac- 
tical purposes (see Ref. |43 for detail benchmarks). In 
RI-V, the bare two-electron integrals are computed as 

{i3\ah) = Y.^^a\^l)V-\v\^h) (A2) 

where 
and 

Here ipp are canonical single-particle spin-orbitals, and 
P^(r) are a set of suitably constructed auxiliary basis 
functions4^ For notational simplicity all orbitals are as- 
sumed to be real. 

In practice, we decompose the matrix in Eq. (jA2[) 
into the product of its square roots, and combine each 
three- index integral with a square root. This gives 



We thank Joachim Paier for making available his 
SOSEX numbers generated using GAUSSIAN, and 



{^J\ab) ^Y.O'iaO';, 



(A5) 
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with 



(A6) 



As discussed in the context of the GW implementation 
in FHI-aims,— the "RI-V" technique can be used to treat 
the screened two-electron Coulomb integrals as well. In 
this case we have 

{ij\W{iio)\ab) = 0^a£;i'{^^)0^, (A7) 



where £ is the coupling-constant averaged dielectric func- 
tions, formally linked to the screened Coulomb matrix by 



-1/2 



(A8) 



In Eq. (jA8|) . W{iuj) is the screened Coulomb interaction 
matrix represented in terms of the auxiliary basis set, 



dvdv'P^{v)W^,{v,v\iu:)P,{v'). (A9) 



For convenience, we introduce a quantity Wiiuj) = 
w^/^Xo(*'^)^^^^^i where xo(*'^) is the independent density 
response function defined in Eq. Using Eqs. ([2]), 

(|A3p . and (|A6p . one can easily obtain the matrix repre- 
sentation of n(iw) in the auxiliary basis 



n (iu) - V '^-^ O^" 



(AlO) 



where and Ca are occupied and unoccupied single- 
particle orbital energies, respectively. Using Eq. (|14p . 
the matrix form of becomes 



dX[l - Xn{iuj)]'^ X. 



(All) 



The A-integration in Eq. (jAlip can be accurately com- 
puted using a Gauss-Legendre quadrature with 5-6 grid 
points. 

Combining Eqs ([XT]) . ((X5|) . and (K7\ . the final expres- 
sion for the RI-SOSEX energy is 



e: 



SOSEX 



° ij,ab L \ M / \ 1^7 



(A12) 



r 



The computational effort for evaluating Eq. (jA12p for- Here is any mean-field potential, which can be non- 
mally scales as 0{N^), where N is the system size. local, as in the case of Hartree-Fock (HF) theory, or local, 

as in the case of Kohn-Sham (KS) theory. 



Appendix B: Derivation of the single excitation 
contribution to the 2nd-order correlation energy 



Suppose the solution of the single-particle Hamiltonian 
is known 



In this section we derive Eq. (|22|) that is presented 
in the main part of this paper - the single excitation 
contribution to the 2nd-order correlation energy ~ from 
Rayleigh-Schrodinger perturbation theory (RSPT). The 
interacting A^-electron system at hand is governed by the 
Hamiltonian 



N 



2^ 



N 

j<k 



1 



rk 



where 1)0x1(1") is a local, multiplicative external potential. 
In RSPT, H is partitioned into a non-interacting mean- 
field Hamiltonian H'^ and an interacting perturbation iJ', 



N N 

= E^°o-) = E 

N 

j<k I 



2^ 



N 

E^r- 



t(r,) 



(Bl) 



then the solution of the non-interacting many-body 
Hamiltonian iJ" follows directly 



i^"l$„) = £;W|a>„). 



The |(E>„) arc single Slater determinants formed from 
of the spin orbitals \p) = j^/'p) determined in Eq. (|Bip . 
These Slater determinants can be distinguished according 
to their excitation level: the ground-state configuration 
|$o), singly excited configurations 1$"), doubly excited 
configurations etc., where i, j, fc, . . . denotes occu- 

pied orbitals and a,h,c, . . . unoccupied ones. Following 
standard perturbation theory, the single-excitation (SE) 
contribution to the 2nd-order correlation energy is given 
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by 



i a ^0 ^ 



EE 



where we have used the fact Eq^^ — E^^^ = a — e^. 



(B2) 



To proceed, the numerator of Eq. (jB2p needs to be 
evaluated. This can most easily be done using second- 
quantization 

^ 1 1 

E |r- -rfcl ^ ^E^P'^I'^*)^^^-^-' 



i<k 



N 



E*r ^ E(pi'^''"i'?)4^.' 

i=j pq 

where p, q, r, s are arbitrary spin-orbitals from Eq. (|Bip . 
c|, and Cq, etc. are the electron creation and annihilation 
operators, and {pq\rs) the two-electron Coulomb inte- 
grals 

',* /'».'\ 



{pq\rs) = j dvdv 



r — r' 



The expectation value of the two-particle Coulomb oper- 
ator between the ground-state configuration <J>o and the 
single excitation $f is given by 

^ occ 

{^Q\2'^{P<lVs)clc\csCr\<^°:) = ^[{ip\ap) - {ip\pa)] 

pqrs p 

= mv^'^l^'a) (B3) 

where v^^ is the HF single-particle potential. 

The expectation value of the mean-field single-particle 
operator , on the other hand, is given by 

($o| J2{p\v'^^\q)clc,\<^>f) = (V'.l^^'^l^a) (B4) 



Combining Eqs. p32|) . p33|) . and pM)) . one gets 



- EE' 



EE 



I Aw,, 



(B5) 

ti — ta 

J o 

where Awia is the matrix element of the difference be- 
tween the HF potential ii^^ and the single-particle mean- 
field potential in question. 

Observing that the ip's are eigenstates of ft." = 
— ^V^ + Vcxt + , and hence all non-diagonal ele- 
ments (^/'i|/i'^|'0a) arc zero, one can alternatively express 
Eq. ^ as 



EE 



KV'.l/IV'a)!^ 



(B6) 



where / is the single-particle HF Hamiltonian, or sim- 
ply Fock operator. Thus Eq. ()22p in the main paper is 
derived. 



For the HF reference state, i.e., when v 



MF 



r,HF 



the tp^s are eigenstates of the Fock operator, and hence 
Eq. (|B2p is zero. For any other reference state, e.g., a KS 
reference state, the ip's are no longer eigenstates of the 
Fock operator, and Eq. (|B2[) is in general not zero. This 
gives rise to a finite SE contribution to the second-order 
correlation energy. 



Appendix C: Derivation of the renormalized single 
excitation (rSE) contribution 

We start with the expression for the second-order 
single-excitation (SE) contribution discussed in Ap- 
pendix |b] 



^SE^V- (*0|i^'|1>.?>(1>.?|i^'|«'0> 



E, 



(0) 



E. 



(0) 



(CI) 



The form of this equation actually already implies that 
the singly excited states |$°) are Slater determinants 
composed of canonical orbitals, namely |^'°) = Det {ipq} 

where = eglV-?), and Ho\<^^} = E^^^l^f) with 



Eia'' = E^"' +<^a-ei- Appendix m shows that Eq. ([Cl 



p(0) 



can be reduced to the simple expression in Eq ()B6P that 
is given in terms of (canonical) single-particle orbitals. 

To set the stage for later discussions, we can also more 
generally express the SE energy in Eq. (jCip in terms of 
non-canonical orbitals {xq}, where ft°|Xp) = J2q ^pq\Xq)j 
and = (XpI" IXg)- In this case, E^^ is given by 



Er = E(<foii?'i$n(*"i(i?r -^o)->?)($?iff'ii>o) 

ij.ab 

= E(^*I/I^'^) [(eI^'^I-Ho)-' 



ia,jb 



iXblflX]) , 



(C2) 



where / is the identity matrix: lia.jb ~ ^ijSat, and 



4°^/ - Ho 



ia.jb 



(C3) 



Now the question arises how to sum up all the higher- 
order SE diagrams shown in Fig. [2]? For canonical or- 
bitals, the corresponding algebraic expression can be eas- 
ily obtained by applying the rules of evaluating Goldstone 
diagrams 
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ia 

+ E 

ijk,c 



faifia 



E 



fai^'^ij fja 



E' 



faiAVtk^Vkjfja 



E 



ij.ab 



E 

ij,ab 



faiA.V^jfjbAvba 



(C4) 



where Avr, 



{^p\f - holipg), and Avia = U 



Eq (|C5p is carried out, we rearrange the expression as 



(V'il/IV'a)- To see how the infinite-order summation in fohows: 



ij.ab 

E 

ijk^abc 
tj,ab 



faiSijSabfjb \ - /ai (AVab^^j - Awij^a;,) /jfc 



E 



/a» (AuifcAwfej^aclJbc + AvbcAVcaSikSkj - AvikAvbcSjkSac ~ Avkj AVcaS^kSbc) f]b 
fai 



(ej - ea)(efc - ec)(ej - Cfc) 



ij.ab 



fai 



^ia^jb ^ia,jb H~ ' ' ' ] /"jb 



where we have introduced the VL matrix, defined as 

AvbaSij - AvijSab 



n 



ia.jb 



£.7 - £6 



(C6) 



Further denoting Aiajb ~ (cj — ^a)SijSab, one observes 
1 



Ci - ea 



[(i-^r'ha,b = [A-'iii-^r'ha,b 

= [iA-^A)-^U^^, (C7) 



and 



{A - fiA)ia,j6 = (f j ea)%'5a(, + AvijSab - AvbaSij 

= /yf^ab — fab^ij , (C8) 

where /y = 6^(5^ + Au,.,, fab = ^aSab + Avab have been 
used. It fohows that 



(C9) 



ij.ab 



We observe that the rSE energy expressed in terms of 
canonical orbitals via Eqs. (|C8[) and (|C9p has the same 
mathematical structure as the second-order SE energy 



(C5) 



r 



expressed in terms of non- canonical orbitals given by 
Eq. (jC2p and (|C3[) . The difference is that now the corre- 
sponding matrix elements in the denominator are evalu- 
ated using the Fock operator /, instead of the KS Hamil- 
tonian operator /i*^. 

To simplify the evaluation of Eq. (jC9|) . one can rotate 
the occupied orbitals and unoccupied orbitals separately, 
such that the Fock matrix becomes diagonal in the occu- 
pied and unoccupied subspaces. This procedure is called 
semi-canonicalization. To be more precise, suppose there 
are transformation matrices O and U which diagonalize 
the fij and fab blocks separately 



^fikOkj = Oijij 

k 

faMcb = l^abh ■ 



(CIO) 



We then have 



J2 OtkKAA - nA)kcMOijUdb = S,jSab{ij - h) (Cll) 



kl.cd 
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3b 



(C14) 



Thus the final expression for rSE has the same form as 
that for SE, only the eigenvalues €^,£0 and the "transi- 
tion amplitude" fia have to be reinterpreted. The actual 
implementation following Eqs. (jClOp . (jClSp . and (|C14|) 
is straightforward. 



FIG. 10: (Color online) The rPT2@PBE binding energy 
of methane dimer in its equilibrium geometry as a function 
of the basis set size. "XZ" and "aXZ" (X=D,T,Q,5) denote 
respectively the Dunning "cc-pVXZ" and "aug-cc-pVXZ" ba- 
sis, whereas "tN" denotes the FHI-aims "tier N" basis, and 
"t3/4" here means tier 4 basis for C and tier 3 basis for H 
(note that a tier 4 basis for H is not available). "t3/4-)-a5Z- 
d" corresponds to the NAO "tier 3/4" plus diffuse functions 
from aug-cc-pV5Z. The ESSE is corrected. 



or equivalently, 

[{A - nA)-']^^^^^ = J2 0..kUac{lk - ~ec)-'Ol^U:, . 

k,c 

(C12) 

Inserting Eq. (IC12P into Eq. (|C9p . one arrives at 




(C13) 



Appendix D: Basis convergence 

Figure [J] shows the convergence behavior of the rPT2 
binding energy of the methane dimer (in its equilib- 
rium geometry) with respect to the FHI-aims NAO "tier 
N" basis as well as Dunning's "cc-pVXZ" and "aug-cc- 
pVXZ" basis. The methane dimer is dominated by the 
dispersion interaction, and the so-called "diffuse func- 
tions" are needed to accurately describe this interaction. 
The difference between the "cc-pVXZ" and "aug-cc- 
pVXZ" results highlight the importance of including "dif- 
fuse functions". For methane dimer the "tier N" series 
exhibits a faster convergence than "cc-pVXZ" whereas 
a slower convergence than "aug-cc-pVXZ" for BSSE- 
corrected binding energies. When adding diffuse func- 
tions from aug-cc-pV5Z to "tier 3/4", (called "t3/4+a5Z- 
d" in Fig. [T0|) results of similar quality as the full aug- 
cc-pV5Z basis are obtained. 
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